void FacetStrain::mask_displace_gradient_filter(
	cv::Mat& padding_lagrange_gradient_xx,
	cv::Mat& padding_lagrange_gradient_yy,
	cv::Mat& padding_lagrange_gradient_xy,
	cv::Mat& padding_euler_gradient_xx,
	cv::Mat& padding_euler_gradient_yy,
	cv::Mat& padding_euler_gradient_xy,
	cv::Mat& padding_mask,
	int halfFilterSize,
	double sigma,
	cv::Mat& result_lagrange_gradient_xx,
	cv::Mat& result_lagrange_gradient_yy,
	cv::Mat& result_lagrange_gradient_xy,
	cv::Mat& result_euler_gradient_xx,
	cv::Mat& result_euler_gradient_yy,
	cv::Mat& result_euler_gradient_xy,
	bool need_lagrange,
	bool need_euler,
	bool crack_mode) {

	double scale2X = -0.5 / (sigma * sigma);
	int filterSize = 2 * halfFilterSize + 1;
	cv::Mat fullFilter = cv::Mat::zeros(filterSize, filterSize, CV_64F);
	cv::Mat unfullFilter = cv::Mat::zeros(filterSize, filterSize, CV_64F);

	// construct full subset filter 
	double distance, weight, weightSum = 0;
	for (int dy = -halfFilterSize; dy <= halfFilterSize; ++dy) {
		for (int dx = -halfFilterSize; dx <= halfFilterSize; ++dx) {
			distance = dy * dy + dx * dx;
			weight = std::exp(scale2X * distance);
			fullFilter.at<double>(dy + halfFilterSize, dx + halfFilterSize) = weight;
			weightSum += weight;
		}
	}
	for (int dy = -halfFilterSize; dy <= halfFilterSize; ++dy) {
		for (int dx = -halfFilterSize; dx <= halfFilterSize; ++dx) {
			fullFilter.at<double>(dy + halfFilterSize, dx + halfFilterSize) /= weightSum;
		}
	}

	cv::Mat filter;
	int padding_height = padding_mask.rows;
	int padding_width = padding_mask.cols;
	cv::Mat subset_mask = cv::Mat::zeros(filterSize, filterSize, CV_8U);

	for (int y = halfFilterSize; y < padding_height - halfFilterSize; ++y) {
		for (int x = halfFilterSize; x < padding_width - halfFilterSize; ++x) {
			if (padding_mask.at<uchar>(y, x) == FILTER_SUBSET_FULL)
				filter = fullFilter;
			else if (padding_mask.at<uchar>(y, x) == FILTER_SUBSET_UNFULL) {
				weightSum = 0;
				cv::Rect subset_rect = cv::Rect(x - halfFilterSize, y - halfFilterSize, filterSize, filterSize);
				if (crack_mode)
					StereoStrainMethod::get_valid_subset(padding_mask, subset_rect, subset_mask);
				else
					subset_mask = padding_mask(subset_rect);
				
				for (int r = 0; r < filterSize; ++r) {
					for (int c = 0; c < filterSize; ++c) {
						if (subset_mask.at<uchar>(r, c) == 0)
							unfullFilter.at<double>(r, c) = 0;
						else {
							int dy = r - halfFilterSize;
							int dx = c - halfFilterSize;
							distance = dy * dy + dx * dx;
							weight = std::exp(scale2X * distance);
							unfullFilter.at<double>(r, c) = weight;
							weightSum += weight;
						}
					}
				}
				for (int r = 0; r < filterSize; ++r)
					for (int c = 0; c < filterSize; ++c)
						unfullFilter.at<double>(r, c) /= weightSum;
				filter = unfullFilter;
			}

			if (!filter.empty()) {
				double sum_lagrange_gradient_xx = 0;
				double sum_lagrange_gradient_yy = 0;
				double sum_lagrange_gradient_xy = 0;
				double sum_euler_gradient_xx = 0;
				double sum_euler_gradient_yy = 0;
				double sum_euler_gradient_xy = 0;
				int yindex, xindex;
				for (int dy = -halfFilterSize; dy <= halfFilterSize; ++dy) {
					for (int dx = -halfFilterSize; dx <= halfFilterSize; ++dx) {
						weight = filter.at<double>(dy + halfFilterSize, dx + halfFilterSize);
						yindex = y + dy;
						xindex = x + dx;
						if (need_lagrange) {
							sum_lagrange_gradient_xx += weight * padding_lagrange_gradient_xx.at<double>(yindex, xindex);
							sum_lagrange_gradient_yy += weight * padding_lagrange_gradient_yy.at<double>(yindex, xindex);
							sum_lagrange_gradient_xy += weight * padding_lagrange_gradient_xy.at<double>(yindex, xindex);
						}
						if (need_euler) {
							sum_euler_gradient_xx += weight * padding_euler_gradient_xx.at<double>(yindex, xindex);
							sum_euler_gradient_yy += weight * padding_euler_gradient_yy.at<double>(yindex, xindex);
							sum_euler_gradient_xy += weight * padding_euler_gradient_xy.at<double>(yindex, xindex);
					